#################################################
# Replication Code
# Coercion and the Credibility of Assurances
# Matthew Cebul, Allan Dafoe, and Nuno Monteiro
#################################################

require(foreign)
require(ggplot2)
require(reshape2)
require("stargazer")
require(plyr)
library(dplyr)
library(sandwich)
library(scales)
require(zeligverse)
require(mediation)
require(car)

set.seed(1991)

# Load Data 
setwd("/Users/mcebul/Dropbox/Co-Authored Projects/Assurances in IR/Submissions/Replication Files")
data <- read.csv("Assurances_JOP_Clean.csv", colClasses=c("Race"="factor"))
attach(data)


########################
# MAIN BODY
########################

##########
# TABLES
##########

#####
# Table 4: Assurance Credibility (Short)
#####

Assurance.1 <- lm(CredAssurance~Power)
Assurance.2 <- lm(CredAssurance~Reputation)
Assurance.3 <- lm(CredAssurance~Power + Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)
Assurance.4 <- lm(CredAssurance~Power + Reputation + Power*Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)

# Robust SEs
Assurance_ses <- lapply(list(Assurance.1, Assurance.2, Assurance.3, Assurance.4), function (x) sqrt(diag(vcovHC(x, type = "HC2"))))

# Construct Table
vars.order <- c("Power", "Reputation", "Power:Reputation", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Edu", "Income", "Male", "MilExp")
stargazer(Assurance.1, Assurance.2, Assurance.3, Assurance.4, title ="Assurance Credibility", 
          se = Assurance_ses,
          dep.var.caption  = "Dependent Variable: Assurance Credibility",
          dep.var.labels   = "If U.S. \\textbf{accepts} demand, will China still contest U.S. military presence?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "Constant", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race"),
          keep=c("Power", "Reputation", "Power:Reputation", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01",
                    "Robust SEs computed via Huber-White sandwhich estimator"), 
          notes.append = FALSE)
# Note: added "Controls" checkmarks + adjusted table spacing manually in LaTeX 


#####
# Table 5: Threat Credibility (Short)
#####

Threat.1 <- lm(CredThreat~Power)
Threat.2 <- lm(CredThreat~Reputation)
Threat.3 <- lm(CredThreat~Power + Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)
Threat.4 <- lm(CredThreat~Power + Reputation + Power*Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)

# Robust SEs
Threat_ses <- lapply(list(Threat.1, Threat.2, Threat.3, Threat.4), function (x) sqrt(diag(vcovHC(x, type = "HC2"))))

# Construct Table
stargazer(Threat.1, Threat.2, Threat.3, Threat.4, title ="Threat Credibility", 
          se = Threat_ses,
          dep.var.caption  = "Dependent Variable: Threat Credibility",
          dep.var.labels   = "If U.S. \\textbf{rejects} demand, will China contest U.S. presence?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "Constant", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race"),
          keep=c("Power", "Reputation", "Power:Reputation", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01",
                    "Robust SEs computed via Huber-White sandwhich estimator"), 
          notes.append = FALSE)
# Note: added "Controls" checkmarks row + adjusted table spacing manually in LaTeX. 


#####
# Table 6: Resolve (Short)
#####

Resolve.1 <- lm(Resolve~Power)
Resolve.2 <- lm(Resolve~Reputation)
Resolve.3 <- lm(Resolve~Power + Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)
Resolve.4 <- lm(Resolve~Power + Reputation + Power*Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)

# Robust SEs
Resolve_ses <- lapply(list(Resolve.1, Resolve.2, Resolve.3, Resolve.4), function (x) sqrt(diag(vcovHC(x, type = "HC2"))))

# Construct Table
stargazer(Resolve.1, Resolve.2, Resolve.3, Resolve.4, title ="Resolve", 
          se = Resolve_ses,
          dep.var.caption  = "Dependent Variable: Resolve",
          dep.var.labels   = "How should U.S. respond to China's demand?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "Constant", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race"),
          keep=c("Power", "Reputation", "Power:Reputation", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01",
                    "Robust SEs computed via Huber-White sandwhich estimator"), 
          notes.append = FALSE)
# Note: added "Controls" checkmarks row + adjusted table spacing manually in LaTeX. 


#####
# Simulating Effect of Rep. Restraint on Resolve
#####

# On pg. 24, we state that a reputation for restraint 
# decreases respondent Resolve to resist Chinese demands by 
# about 13%, holding all other parameters at their means. This 
# result was obtained via simulation, using Zelig:

# Full model
z.model <- zelig(Resolve~Power + Reputation + GOPPresident + Age
                 + Race + Dem + Rep + PolKnow + Edu + Income + Male
                 + MilExp, model = 'ls', data = data)

# Set Reputation Values
x.replow <- setx(z.model, Reputation = 0, Race = 1)
x.rephigh <- setx(z.model, Reputation = 1, Race = 1)

# Simulate
s.outrep <- sim(z.model, x = x.replow, x1 = x.rephigh, num = 1000000)

# Calculate Percent Change
summary(s.outrep)
# E(Resolve | Rep. Restraint = 0) = 1.87
# E(Resolve | Rep. Restraint = 1) = 1.62
# (1.87 - 1.62)/1.87 = approximately 13% decrease


##########
# FIGURES
##########

# Figures 2 and 3 re-display data from the regression
# tables -- replicating them is an exercise in ggplot 
# figure construction. For replication code for
# the mediation analysis, skip to Figure 4.

#####
# Figure 2, Top: Assurance Credibility
#####

# Temporary data.frames:
AssuranceF1 <- data.frame(Variable = rownames(summary(Assurance.1)$coef),
                          Coefficient = summary(Assurance.1)$coef[, 1],
                          SE = summary(Assurance.1)$coef[, 2],
                          modelName = "Univariate, N = 1,028")
AssuranceF2 <- data.frame(Variable = rownames(summary(Assurance.2)$coef),
                          Coefficient = summary(Assurance.2)$coef[, 1],
                          SE = summary(Assurance.2)$coef[, 2],
                          modelName = "Univariate, N = 1,028")
AssuranceF3 <- data.frame(Variable = rownames(summary(Assurance.3)$coef),
                          Coefficient = summary(Assurance.3)$coef[, 1],
                          SE = summary(Assurance.3)$coef[, 2],
                          modelName = "Controls, N = 970")
AssuranceF4 <- data.frame(Variable = rownames(summary(Assurance.4)$coef),
                              Coefficient = summary(Assurance.4)$coef[, 1],
                              SE = summary(Assurance.4)$coef[, 2],
                              modelName = "Interacted, N = 970")

# Combine them
AssuranceFull <- data.frame(rbind(AssuranceF1, AssuranceF2, AssuranceF3, AssuranceF4))

# Add robust SEs
credassurance_se_frame <- do.call(rbind, lapply(Assurance_ses, data.frame, stringsAsFactors=FALSE))
AssuranceFull$SE <- credassurance_se_frame$X..i..

#Drop intercepts
AssuranceFull <- AssuranceFull[AssuranceFull$Variable != "(Intercept)",] 

# Change Labels + Order
AssuranceFull$Variable <- factor(AssuranceFull$Variable, levels = c("Power", "Reputation", "GOPPresident", "Power:Reputation", "Dem", "Rep", "PolKnow", "Edu", "Income", "Male", "MilExp", "Age", "Race2", "Race3", "Race4", "Race7"))
AssuranceFull$Variable <- revalue(AssuranceFull$Variable, c("Power" = "China Power", "Reputation" = "Reputation", "GOPPresident"="GOP President", "Power:Reputation" = "Power*Reputation"))
colnames(AssuranceFull)[4] <- "ModelType"

# Drop Controls
AssuranceKey <- AssuranceFull[AssuranceFull$Variable == "China Power" | AssuranceFull$Variable == "Reputation" | AssuranceFull$Variable == "Power*Reputation",]
AssuranceKey$Variable <- factor(AssuranceKey$Variable)

# Confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
p1 <- ggplot(AssuranceKey, aes(colour = ModelType)) + scale_colour_grey(start = 0, end = .75)
p1 <- p1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
p1 <- p1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
p1 <- p1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
p1 <- p1 + theme_bw() + scale_x_discrete(limits = rev(levels(AssuranceKey$Variable))) + theme(text = element_text(size=14))
p1 <- p1 + coord_flip() + theme(aspect.ratio = 1)
p1 <- p1 + ggtitle("Assurance Credibility") + theme(plot.title = element_text(hjust = 0.5))
print(p1) 

#####
# Figure 2, Bottom: Threat Credibility
#####

# Temporary data.frames
ThreatF1 <- data.frame(Variable = rownames(summary(Threat.1)$coef),
                          Coefficient = summary(Threat.1)$coef[, 1],
                          SE = summary(Threat.1)$coef[, 2],
                          modelName = "Univariate, N = 1,028")
ThreatF2 <- data.frame(Variable = rownames(summary(Threat.2)$coef),
                          Coefficient = summary(Threat.2)$coef[, 1],
                          SE = summary(Threat.2)$coef[, 2],
                          modelName = "Univariate, N = 1,028")
ThreatF3 <- data.frame(Variable = rownames(summary(Threat.3)$coef),
                          Coefficient = summary(Threat.3)$coef[, 1],
                          SE = summary(Threat.3)$coef[, 2],
                          modelName = "Controls, N = 970")
ThreatF4 <- data.frame(Variable = rownames(summary(Threat.4)$coef),
                              Coefficient = summary(Threat.4)$coef[, 1],
                              SE = summary(Threat.4)$coef[, 2],
                              modelName = "Interacted, N = 970")

# Combine them
ThreatFull <- data.frame(rbind(ThreatF1, ThreatF2, ThreatF3, ThreatF4))

# Add robust SEs
credthreat_se_frame <- do.call(rbind, lapply(Threat_ses, data.frame, stringsAsFactors=FALSE))
ThreatFull$SE <- credthreat_se_frame$X..i..

#Drop intercepts
ThreatFull <- ThreatFull[ThreatFull$Variable != "(Intercept)",] 

# Change Labels and Order
ThreatFull$Variable <- factor(ThreatFull$Variable, levels = c("Power", "Reputation", "GOPPresident", "Power:Reputation", "Dem", "Rep", "PolKnow", "Edu", "Income", "Male", "MilExp", "Age", "Race2", "Race3", "Race4", "Race7"))
ThreatFull$Variable <- revalue(ThreatFull$Variable, c("Power" = "China Power", "Reputation" = "Reputation", "GOPPresident"="GOP President", "Power:Reputation" = "Power*Reputation"))
colnames(ThreatFull)[4] <- "ModelType"

# Drop Controls
ThreatKey <- ThreatFull[ThreatFull$Variable == "China Power" | ThreatFull$Variable == "Reputation" | ThreatFull$Variable == "Power*Reputation",]
ThreatKey$Variable <- factor(ThreatKey$Variable)

# Confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
p2 <- ggplot(ThreatKey, aes(colour = ModelType)) + scale_colour_grey(start = 0, end = .75)
p2 <- p2 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
p2 <- p2 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
p2 <- p2 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
p2 <- p2 + theme_bw() + scale_x_discrete(limits = rev(levels(ThreatKey$Variable))) + theme(text = element_text(size=14))
p2 <- p2 + coord_flip() + theme(aspect.ratio = 1)
p2 <- p2 + ggtitle("Threat Credibility") + theme(plot.title = element_text(hjust = 0.5))
print(p2) 

#####
# Figure 3: Resolve
#####

# Temporary data.frames:
ResolveF1 <- data.frame(Variable = rownames(summary(Resolve.1)$coef),
                          Coefficient = summary(Resolve.1)$coef[, 1],
                          SE = summary(Resolve.1)$coef[, 2],
                          modelName = "Univariate, N = 1,028")
ResolveF2 <- data.frame(Variable = rownames(summary(Resolve.2)$coef),
                          Coefficient = summary(Resolve.2)$coef[, 1],
                          SE = summary(Resolve.2)$coef[, 2],
                          modelName = "Univariate, N = 1,028")
ResolveF3 <- data.frame(Variable = rownames(summary(Resolve.3)$coef),
                        Coefficient = summary(Resolve.3)$coef[, 1],
                        SE = summary(Resolve.3)$coef[, 2],
                        modelName = "Controls, N = 970")
ResolveF4 <- data.frame(Variable = rownames(summary(Resolve.4)$coef),
                            Coefficient = summary(Resolve.4)$coef[, 1],
                            SE = summary(Resolve.4)$coef[, 2],
                            modelName = "Interacted, N = 970")

# Combine them
ResolveFull<- data.frame(rbind(ResolveF1, ResolveF2, ResolveF3, ResolveF4))

# Add robust SEs
resolve_se_frame <- do.call(rbind, lapply(Resolve_ses, data.frame, stringsAsFactors=FALSE))
ResolveFull$SE <- resolve_se_frame$X..i..

#Drop intercepts
ResolveFull <- ResolveFull[ResolveFull$Variable != "(Intercept)",] 

# Change Labels and Order
ResolveFull$Variable <- factor(ResolveFull$Variable, levels = c("Power", "Reputation", "GOPPresident", "Power:Reputation", "Dem", "Rep", "PolKnow", "Edu", "Income", "Male", "MilExp", "Age", "Race2", "Race3", "Race4", "Race7"))
ResolveFull$Variable <- revalue(ResolveFull$Variable, c("Power" = "China Power", "Reputation" = "Reputation", "GOPPresident"="GOP President", "Power:Reputation" = "Power*Reputation"))
colnames(ResolveFull)[4] <- "ModelType"

# Drop Controls
ResolveKey <- ResolveFull[ResolveFull$Variable == "China Power" | ResolveFull$Variable == "Reputation" | ResolveFull$Variable == "Power*Reputation",]
ResolveKey$Variable <- factor(ResolveKey$Variable)

# Confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
p3 <- ggplot(ResolveKey, aes(colour = ModelType)) + scale_colour_grey(start = 0, end = .75)
p3 <- p3 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
p3 <- p3 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
p3 <- p3 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
p3 <- p3 + theme_bw() + scale_x_discrete(limits = rev(levels(ResolveKey$Variable))) + theme(text = element_text(size=14))
p3 <- p3 + coord_flip() + theme(aspect.ratio = 1)
p3 <- p3 + ggtitle("Resolve") + theme(plot.title = element_text(hjust = 0.5))
print(p3) 

#####
# Figure 4: Mediation Analysis
#####

# Threat and assurance credibility are correlated mediators. Thus,
# we employ Imai and Yamamoto's 2013 procedure for mediation
# analysis with multiple, correlated mediators:

# Power, Primary Mediator = CredAssurance
powermed <- multimed(outcome = "Resolve", med.main = "CredAssurance", 
                     med.alt = "CredThreat", treat = "Power", 
                     covariates = c("Reputation", "GOPPresident", "Age", "Race", "Dem", "Rep", "PolKnow", "Edu", "Income", "Male", "MilExp"),
                     sims = 10000, conf.level = 0.95, data = data)

# Power, Primary Mediator = CredThreat
powermed2 <- multimed(outcome = "Resolve", med.main = "CredThreat", 
                      med.alt = "CredAssurance", treat = "Power", 
                      covariates = c("Reputation", "GOPPresident", "Age", "Race", "Dem", "Rep", "PolKnow", "Edu", "Income", "Male", "MilExp"),
                      sims = 10000, conf.level = 0.95, data = data)

# Reputation, Primary Mediator = CredAssurance
assurancemed <- multimed(outcome = "Resolve", med.main = "CredAssurance", 
                         med.alt = "CredThreat", treat = "Reputation", 
                         covariates = c("Power", "GOPPresident", "Age", "Race", "Dem", "Rep", "PolKnow", "Edu", "Income", "Male", "MilExp"),
                         sims = 10000, conf.level = 0.95, data = data)

# Reputation, Primary Mediator = CredThreat
assurancemed2 <- multimed(outcome = "Resolve", med.main = "CredThreat", 
                          med.alt = "CredAssurance", treat = "Reputation", 
                          covariates = c("Power", "GOPPresident", "Age", "Race", "Dem", "Rep", "PolKnow", "Edu", "Income", "Male", "MilExp"),
                          sims = 10000, conf.level = 0.95, data = data)

# Plot them 
par(mfrow = c(2,2))
par(mar=c(5,6.8,4,2))
PowerAssurance <- plot(powermed, type = "point", tgroup = "average", 
                                xlim = c(-.4, .15), 
                                main = expression(Power %->% CredAssurance %->% Resolve), pch = 19, 
                                xlab = "Average Effects", ylab = c("Total Effect", "Direct Effect", "CredAssurance Effect"))
PowerThreat <- plot(powermed2, type = "point", tgroup = "average", 
                             xlim = c(-.4, .15), 
                             pch = 19, main = expression(Power %->% CredThreat %->% Resolve),
                             xlab = "Average Effects", ylab = c("Total Effect", "Direct Effect", "CredThreat Effect"))
RepAssurance <- plot(assurancemed, type = "point", tgroup = "average", 
                              xlim = c(-.4, .15), 
                              pch = 19, main = expression(Rep. %->% CredAssurance %->% Resolve),
                              xlab = "Average Effects", ylab = c("Total Effect", "Direct Effect", "CredAssurance Effect"))
RepThreat <- plot(assurancemed2, type = "point", tgroup = "average", 
                               xlim = c(-.4, .15), 
                               pch = 19, main = expression(Rep. %->% CredThreat %->% Resolve),
                               xlab = "Average Effects", ylab = c("Total Effect", "Direct Effect", "CredThreat Effect"))


########################
# APPENDIX
########################


#####
# Section 1: Full Regression Tables 
#####

# Table 1: Assurance Credibility (Full)
vars.order <- c("Power", "Reputation", "Power:Reputation", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Edu", "Income", "Male", "MilExp")
stargazer(Assurance.1, Assurance.2, Assurance.3, Assurance.4, title="Assurance Credibility", 
          se = Assurance_ses,
          dep.var.caption  = "Dependent Variable: Assurance Credibility",
          dep.var.labels   = "If U.S. \\textbf{accepts} demand, will China still contest U.S. military presence?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01",
                    "Robust SEs computed via Huber-White sandwhich estimator"), 
          notes.append = FALSE)

# Table 2: Threat Credibility (Full)
stargazer(Threat.1, Threat.2, Threat.3, Threat.4, title="Threat Credibility", 
          se = Threat_ses,
          dep.var.caption  = "Dependent Variable: Threat Credibility",
          dep.var.labels   = "If U.S. \\textbf{rejects} demand, will China militarize?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01",
                    "Robust SEs computed via Huber-White sandwhich estimator"), 
          notes.append = FALSE)

# Table 3: Resolve (Full)
stargazer(Resolve.1, Resolve.2, Resolve.3, Resolve.4, title ="Resolve", 
          se = Resolve_ses,
          dep.var.caption  = "Dependent Variable: Resolve",
          dep.var.labels   = "How should U.S. respond to China's demand?",
          float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01",
                    "Robust SEs computed via Huber-White sandwhich estimator"), 
          notes.append = FALSE)


#####
# Section 2: Balance Tests
#####

# Making Race Dummies
White <- ifelse(Race==1, 1, 0)
AfAm <- ifelse(Race == 2, 1, 0)
Asian <- ifelse(Race==3, 1, 0)
Hispanic <- ifelse(Race==4, 1, 0)
OtherRace <- ifelse(Race==7,1,0)

# Balance Table, Power
balance.pow <- as.data.frame(rbind(
  summary(lm((Age/10)~Power))$coef[2,1:2],
  summary(lm(White~Power))$coef[2,1:2],
  summary(lm(AfAm~Power))$coef[2,1:2],
  summary(lm(Asian~Power))$coef[2,1:2],
  summary(lm(Hispanic~Power))$coef[2,1:2],
  summary(lm(OtherRace~Power))$coef[2,1:2],
  summary(lm(Dem~Power))$coef[2,1:2],
  summary(lm(Rep~Power))$coef[2,1:2],
  summary(lm(PolKnow~Power))$coef[2,1:2],
  summary(lm(Edu~Power))$coef[2,1:2],
  summary(lm(Income/2~Power))$coef[2,1:2],
  summary(lm(Male~Power))$coef[2,1:2],
  summary(lm(MilExp~Power))$coef[2,1:2]
))

# Balance Table, Reputation
balance.rep <- as.data.frame(rbind(
  summary(lm((Age/10)~Reputation))$coef[2,1:2],
  summary(lm(White~Reputation))$coef[2,1:2],
  summary(lm(AfAm~Reputation))$coef[2,1:2],
  summary(lm(Asian~Reputation))$coef[2,1:2],
  summary(lm(Hispanic~Reputation))$coef[2,1:2],
  summary(lm(OtherRace~Reputation))$coef[2,1:2],
  summary(lm(Dem~Reputation))$coef[2,1:2],
  summary(lm(Rep~Reputation))$coef[2,1:2],
  summary(lm(PolKnow~Reputation))$coef[2,1:2],
  summary(lm(Edu~Reputation))$coef[2,1:2],
  summary(lm(Income/2~Reputation))$coef[2,1:2],
  summary(lm(Male~Reputation))$coef[2,1:2],
  summary(lm(MilExp~Reputation))$coef[2,1:2]
))

# Adjust labels 
balance.pow$name <- factor(c("Age/10", "Caucasian (0-1)", "African American (0-1)", "Asian (0-1)", "Hispanic (0-1)", "OtherRace (0-1)", "Dem (0-1)", "Rep (0-1)", "PolKnow (0-4)","Education (0-6)", "Income (0-5)","Male (0-1)", "MilExperience (0-1)"),
                           levels = c("Age/10", "Caucasian (0-1)", "African American (0-1)", "Asian (0-1)", "Hispanic (0-1)", "OtherRace (0-1)", "Dem (0-1)", "Rep (0-1)", "PolKnow (0-4)","Education (0-6)", "Income (0-5)","Male (0-1)", "MilExperience (0-1)"))
balance.rep$name <- factor(c("Age/10", "Caucasian (0-1)", "African American (0-1)", "Asian (0-1)", "Hispanic (0-1)", "OtherRace (0-1)", "Dem (0-1)", "Rep (0-1)", "PolKnow (0-4)","Education (0-6)", "Income (0-5)","Male (0-1)", "MilExperience (0-1)"),
                           levels = c("Age/10", "Caucasian (0-1)", "African American (0-1)", "Asian (0-1)", "Hispanic (0-1)", "OtherRace (0-1)", "Dem (0-1)", "Rep (0-1)", "PolKnow (0-4)","Education (0-6)", "Income (0-5)","Male (0-1)", "MilExperience (0-1)"))
colnames(balance.pow) <- c("mean", "se", "name")
colnames(balance.rep) <- c("mean", "se", "name")

# Plot Them
bp <- ggplot(balance.pow)
bp <- bp + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
bp <- bp + geom_pointrange(aes(x = name, y = mean, ymin = mean - se*interval2,
                                 ymax = mean + se*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
bp <- bp + theme_bw() + scale_x_discrete(limits = rev(levels(balance.pow$name)))
bp <- bp + coord_flip(ylim=c(-0.25,0.4)) + labs(y = "Difference in Means", x= "")
bp <- bp + ggtitle("Balance Plot, Power") + theme(plot.title = element_text(hjust = 0.5))
print(bp) 

br <- ggplot(balance.rep)
br <- br + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
br <- br + geom_pointrange(aes(x = name, y = mean, ymin = mean - se*interval2,
                                 ymax = mean + se*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
br <- br + theme_bw() + scale_x_discrete(limits = rev(levels(balance.rep$name)))
br <- br + coord_flip(ylim=c(-0.4,0.3)) + labs(y = "Difference in Means", x= "")
br <- br + ggtitle("Balance Plot, Reputation") + theme(plot.title = element_text(hjust = 0.5))
print(br) 


#####
# Section 3: Bonferroni Corrections 
#####

# Define Bonferroni Function
bonferroni <- function(x) (p.adjust(x, method = "bonferroni", n=4))

# Table 4: Assurance Credibility, Bonferroni
stargazer(Assurance.1, Assurance.2, Assurance.3, Assurance.4, title="Assurance Credibility, Bonferroni", 
          se = Assurance_ses,
          report=('vc*p'),
          apply.p = bonferroni,
          dep.var.caption  = "Dependent Variable: Assurance Credibility",
          dep.var.labels   = "If U.S. \\textbf{accepts} demand, will China still contest U.S. military presence?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "Constant", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race"),
          keep=c("Power", "Reputation", "Power:Reputation", "Constant"),
           notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01"), 
          notes.append = FALSE)
# Note: added "Controls" checkmarks row + adjusted table spacing manually in LaTeX. 

# Table 5: Threat Credibility, Bonferroni
stargazer(Threat.1, Threat.2, Threat.3, Threat.4, title="Threat Credibility, Bonferroni", 
          se = Threat_ses,
          report=('vc*p'),
          apply.p = bonferroni,
          dep.var.caption  = "Dependent Variable: Threat Credibility",
          dep.var.labels   = "If U.S. \\textbf{rejects} demand, will China militarize?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "Constant", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race"),
          keep=c("Power", "Reputation", "Power:Reputation", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01"), 
          notes.append = FALSE)
# Note: added "Controls" checkmarks row + adjusted table spacing manually in LaTeX. 

# Table 6: Resolve, Bonferroni
stargazer(Resolve.1, Resolve.2, Resolve.3, Resolve.4, title="Resolve, Bonferroni", 
          se = Resolve_ses,
          report=('vc*p'),
          apply.p = bonferroni,
          dep.var.caption  = "Dependent Variable: Resolve",
          dep.var.labels   = "How should U.S. respond to China's demand?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "Constant", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race"),
          keep=c("Power", "Reputation", "Power:Reputation", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01"), 
          notes.append = FALSE)
# Note: added "Controls" checkmarks row + adjusted table spacing manually in LaTeX. 


#############
# Section 4: Elite Subset 
#############

data.elite <- data[data$Age >= 40,] #subset to respondents older than 40 years
data.elite <- data.elite[data.elite$Edu >= 4,] #subset to respondents with a college degree
length(data.elite$Resolve) #296 respondents remaining

detach(data)
attach(data.elite)

###
# Table 7: Assurance Credibility, Elites
###

EAssurance.1 <- lm(CredAssurance~Power)
EAssurance.2 <- lm(CredAssurance~Reputation)
EAssurance.3 <- lm(CredAssurance~Power + Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)
EAssurance.4 <- lm(CredAssurance~Power + Reputation + Power*Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)

# Robust SEs
EAssurance_ses <- lapply(list(EAssurance.1, EAssurance.2, EAssurance.3, EAssurance.4), function (x) sqrt(diag(vcovHC(x, type = "HC2"))))

# Construct Table
stargazer(EAssurance.1, EAssurance.2, EAssurance.3, EAssurance.4, title ="Assurance Credibility, Quasi-Elites", 
          se = EAssurance_ses,
          dep.var.caption  = "Dependent Variable: Assurance Credibility",
          dep.var.labels   = "If U.S. \\textbf{accepts} demand, will China still contest U.S. military presence?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "Constant", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race"),
          keep=c("Power", "Reputation", "Power:Reputation", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01",
                    "Robust SEs computed via Huber-White sandwhich estimator"), 
          notes.append = FALSE)

###
# Table 8: Threat Credibility, Elites
### 

EThreat.1 <- lm(CredThreat~Power)
EThreat.2 <- lm(CredThreat~Reputation)
EThreat.3 <- lm(CredThreat~Power + Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)
EThreat.4 <- lm(CredThreat~Power + Reputation + Power*Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)

# Robust SEs
EThreat_ses <- lapply(list(EThreat.1, EThreat.2, EThreat.3, EThreat.4), function (x) sqrt(diag(vcovHC(x, type = "HC2"))))

# Construct Table
stargazer(EThreat.1, EThreat.2, EThreat.3, EThreat.4, title ="Threat Credibility, Quasi-Elites", 
          se = EThreat_ses,
          dep.var.caption  = "Dependent Variable: Threat Credibility",
          dep.var.labels   = "If U.S. \\textbf{rejects} demand, will China contest U.S. presence?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "Constant", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race"),
          keep=c("Power", "Reputation", "Power:Reputation", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01",
                    "Robust SEs computed via Huber-White sandwhich estimator"), 
          notes.append = FALSE)

###
# Table 9: Resolve, Quasi-Elites
###

EResolve.1 <- lm(Resolve~Power)
EResolve.2 <- lm(Resolve~Reputation)
EResolve.3 <- lm(Resolve~Power + Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)
EResolve.4 <-  lm(Resolve~Power + Reputation + Power*Reputation + GOPPresident + Age + Race + Dem + Rep + PolKnow + Edu + Income + Male + MilExp)

# Robust SEs
EResolve_ses <- lapply(list(EResolve.1, EResolve.2, EResolve.3, EResolve.4), function (x) sqrt(diag(vcovHC(x, type = "HC2"))))

# Construct Table
stargazer(EResolve.1, EResolve.2, EResolve.3, EResolve.4, title ="Resolve, Quasi-Elites", 
          se = EResolve_ses,
          dep.var.caption  = "Dependent Variable: Resolve",
          dep.var.labels   = "How should U.S. respond to China's demand?",
          align=TRUE, float = TRUE, column.sep.width = "0pt", 
          star.char = c("\\dagger", "*", "**"),
          order=paste0("^", vars.order , "$"),
          covariate.labels = c("China Power", "China Reputation", "Power*Reputation", "Constant", "GOPPresident", "Age", "Dem", "Rep", "PolKnow", "Education", "Income", "Male", "MilExp", "African-American", "Asian", "Hispanic", "Other Race"),
          keep=c("Power", "Reputation", "Power:Reputation", "Constant"),
          notes = c("$^{\\dagger}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01",
                    "Robust SEs computed via Huber-White sandwhich estimator"), 
          notes.append = FALSE)


###################################################
###################################################
